Relative dispersion in fully developed turbulence: 
The Richardson's Law and Intermittency Corrections 
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Relative dispersion in fully developed turbulence is investigated by means of direct numerical 
simulations. Lagrangian statistics is found to be compatible with Richardson description although 
small systematic deviations are found. The value of the Richardson constant is estimated as C2 — 
0.55, in a close agreement with recent experimental findings [S. Ott and J. Mann J. Fluid Mech. 422, 
207 (2000)]. By means of exit-time statistics it is shown that the deviations from Richardson's law 
are a consequence of Eulerian intermittency. The measured Lagrangian scaling exponents require a 
set of Eulerian structure function exponents £ p which are remarkably close to standard ones known 
for fully developed turbulence. 



The statistics of two particle dispersion is historically 
one of the first issues which has been quantitatively ad- 
dressed in the study of fully developed turbulence. This 
was done by Richardson, in a pioneering work on the 
properties of dispersion in the atmosphere in 1926 El, 15 
years before the theoretical development by Kolmogorov 
and Obukhov H. Despite this fact, there are still rela- 
tively few experimental studies on turbulent Lagrangian 
dispersion. This is essentially due to the difficulties to ob- 
tain Lagrangian trajectories in fully developed turbulent 
flow. The first studies where done in geophysical flows 
(sec [|| for a review) in which Lagrangian tracers are more 
easily followed. Recently, the problem was approached in 
laboratory experiments j|||] but the results are still not 
conclusive. Moreover, most of the numerical studies of 
relative dispersion rely on kinematic simulations in syn- 
thetic flows |5jj6J. Direct numerical simulation have been 
done mostly for two-dimensional turbulence |7]||. 

The scope of this Letter is to contribute to the under- 
standing of relative dispersion by means of direct numer- 
ical simulations of three dimensional turbulence. In what 
follows we show the qualitative validity of the Richard- 
son's description, and discuss its limitations as posed by 
Lagrangian intermittency, whose properties will be inves- 
tigated in detail. 

Richardson's original description of relative disper- 
sion is based on a diffusion equation for the probability 
density function of pair separation p(r, t) which in the 
isotropic case can be written as 



dp(r, t) 
dt 
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dp(r, t) 
dr 



(1) 



The turbulent eddy diffusivity was empirically estab- 
lished by Richardson to follow the "four-thirds law" 
K(r) oc r 4 / 3 . This law is a direct consequence of the 
small-scale velocity statistics, as was first recognized by 
Obukhov P). Thus, for r within the inertial range, the 
dimensional analysis gives 



K{r) = fc e 1/3 r 4 / 3 , 



(2) 



where e is the mean energy dissipation and kg a di- 
mensionless constant. We remark that (0) does not 
imply that a finite energy flux is necessary for parti- 
cle dispersion. Indeed, particle separation is observed 
also in pseudo-turbulent synthetic Gaussian velocity field 
HIJrJ. Using (|), the solution of (@) for ^-distribution 
initial condition has the well known form 



p{r,t) 
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where A = 2187/22407T 3 / 2 is a normalizing factor. The 
most important feature of the Richardson distribution 
(0) is non-Gaussianity with a very pronounced peak at 
the origin and rather fat tails. In the past, alternative 
distributions have been proposed ET 12 1. In particular 
Batchelor ElJ suggested a Gaussian distribution as a con- 
sequence of a diffusivity which depends only on averaged 
quantities. Because the available data is scarce, there is 
still no general consensus on the real form of pair sepa- 
ration pdf. Recent experimental works |0,0| are in favor 
of |). 

The possibility to describe the dispersion process by 
means of a diffusion equation is based on two physical 
assumptions. The first is that the velocity field is short 
correlated in time. Indeed, in the limit of velocity field 
(^-correlated in time (the so-called Kraichnan model of 
turbulence) the diffusion equation of the type of Eq.(Q) 
becomes exact jl3|,[14|. The effects of finite correlation 
time have been recently discussed in p fD^Jlql . 

The second assumption, which is one of the points dis- 
cussed in this Letter, is that the dispersion process is 
self-similar in time, i.e. the scaling exponents of the mo- 
ments of the separation 



R 2n {t) = (r 2n (t)) = C 2n e n t c 



(4) 



have the values a^n = 3n/2, as following from dimen- 
sional analysis. If this is the case, a single number, such 
as the Richardson constant Ci is sufficient to parameter- 
ize turbulent dispersion. There is still a large uncertainty 
on the value of C2, ranging from 0(1O -2 ) — O(IC)™ 1 ) for 
kinematic simulations MM to 0(1) or more in the case 
of closure predictions B. A recent experimental inves- 
tigation gives the value Ci = 0.5 [01. The hypothesis 
of self-similarity is reasonable with a self-affine Eulerian 
velocity, such as in the case of two-dimensional inverse 
cascade turbulence [pj . A recent analysis of a kinematic 
model with synthetic velocity field has shown that La- 
grangian self-similarity can be broken in presence of Eu- 
lerian intermittency. In this case the exponents a n have 
been found in agreement with the prediction of a multi- 
fractal approach for Lagrangian statistics. In particular, 
the second moment of relative dispersion is not affected 
by intermittency, i.e. 011 — 3 fllOfl , essentially because it 
is proportional to e . We remind that Lagrangian inter- 
mittency has been observed also in the case of the so- 
called strong anomalous diffusion |17J . Although in that 
case the mechanism leading to intermittency is different 
(there is no scaling invariant flow), the implication for 
Lagrangian description is identical, i.e. the process can 
not be described by a Fokker-Planck equation like (|J). 

We now turn to our numerical procedure. The tur- 
bulent velocity field is generated by direct integration 
of Navier-Stokes equation in a periodic box of size L — 
27r. The integration is done on a parallel computer by 
means of a pseudo-spectral code at resolution 256 3 with 
Re\ ~ 200. Energy is injected into the flow by keep- 
ing the total energy in each of the two first wavenum- 
ber shells constant J18| and is removed by a second-order 
hyperviscous dissipation. Time integration is performed 
with a second-order Runge-Kutta scheme. In Figure [y 
we plot the energy spectrum which shows a well devel- 
oped Kolmogorov power-law scaling. A small "bump" at 
k > 20 is the signature of a bottleneck effect jl9| as a 
consequence of hyperviscosity. In the inset of Figure [l] 
we show the third order longitudinal structure function 
Ss(x) — (5u(x) 3 ) compensated with the theoretical pre- 
diction S 3 (x) = -4/5ex ||o). 

Passive tracer trajectories are obtained by integrating 
x(£) = u(x(£),£) with the velocity at particle positions 
obtained by linear interpolation from the nearest grid 
points. The reported results are obtained averaging over 
a total number of about 3 x 10 5 particle pairs starting 
from initial separation R(0) = L/256 and over 7 large 
scale eddy turnover times. 

In Figure @ we plot the second moment of relative dis- 
persion R 2 (t). The Richardson t 3 law fl) is clearly ob- 
servable although systematic deviations are detectable, 
in particular in the compensated plot. These deviations, 
observed also in kinematic simulations JlOJ ] and in two- 
dimensional turbulence [pi, are due to finite size effects. 
Consider a series of pair dispersion experiments, in each 



of which a couple of particles is released at time t = at 
initial separation i?(0). At a fixed time t one performs 
an average over all realizations and computes R 2 (t) . For 
t small R 2 (t) is dominated by the initial distance, so 
that the i? 2 (i)-curve flattens. For large times some pairs 
might have reached a separation larger than the integral 
scale and thus show normal (not Richardson) diffusion, 
so that the R 2 (i)-dependence flattens again. Under these 
conditions, a precise determination of the exponents and 
coefficients in (0), in particular the Richardson constant 
C2, is very difficult. 

The distribution of relative separations is plotted in 
Figure || for three different times. The form of the pdf is 
very close to the Richardson prediction (||) and excludes 
other distributions. Our result is the first direct numer- 
ical evidence of the substantial validity of Richardson's 
equation and gives support to recent experimental find- 
ings [fl . A closer inspection of Figure reveals however 
that the self-similar evolution predicted by (pi) is not ex- 
act. Again, the deviations from the distribution (0) are 
mostly due to finite Reynolds effects: because of the large 
tails, a large fraction of particles exits the inertial range 
after a very short time. 

To overcome these difficulties in Lagrangian statistics, 
an alternative approach based on doubling time (or exit 
time) statistics has been recently proposed Jl(],|2l| . Given 
a set of thresholds R n — p n R(Q) within the inertial range, 
one computes the "doubling time" T p (R n ) defined as the 
time it takes for the particle pair separation to grow from 
threshold R n to the next one R n +i- Averages are then 
performed over many dispersion experiments, i.e., par- 
ticle pairs. The outstanding advantage of averaging at 
fixed scale separation, as opposite to a fixed time, is 
that it removes crossover effects since all sampled par- 
ticle pairs belong to the same scales. In the simulations 
presented here, the value p — 1.2 is used. 

Let us first show how doubling time analysis can be 
used for estimating the Richardson constant C%. Ne- 
glecting intermittency, the mean doubling time can be 
obtained from the first-passage problem for the Richard- 
son diffusion equation ( jj) as || 



(T P (R)) 



p 2 ' 3 - 1 
2fc e 1/3 P 2/3 



R 2 ' 3 . 



(5) 



From (0) and (0) one has C2 = ^^k^. Comparison with 
(§) gives 



C 2 = 



143 (p 



2/3 



R 2 



SI 



<W 



(6) 



In the inset of Figure we plot expression (pi) which gives 
directly the value of C 2 . Although the compensation is 
not perfect, it is possible to estimate the Richardson con- 
stant with much better accuracy than from the direct 
analysis of Figure ||. The resulting value, C2 = 0.55±0.1, 



is remarkably close to the recent experimental finding |jj . 
The non perfect compensation is the consequence of in- 
termittency. 

Let us now discuss the issue of intermittency in more 
detail and concentrate on the behavior of the moments 
of inverse doubling times, ((l/T p (R)) p ). We expect for 
doubling time statistics a power-law behavior 



T P {R) 



}~rP* 



(7) 



with exponents p connected to the exponents a n in 
(jfy. Negative moments of doubling time are dominated 
by pairs which separate fast; this corresponds to posi- 
tive moments of relative separation. Kolmogorov scaling, 
based on the dimensional analysis, gives ((1/T P (R)) P ) ~ 
£ p/3ft-2 P /3 g0 foafc p p = —2p/3. Intermittency can be 
taken into account by using the simple dimensional es- 
timate for the doubling time, T(R) ~ R/6u(R) which 
gives 



P P = ( P - p, 



(8) 



where Cp are the scaling exponents of the longitudinal 
structure functions. As a consequence of the Kolmogorov 
"4/5" law, £3 = 1 pTJ and the doubling time exponent 
not affected by intermittency is 0% = —2 (again, the 
quantity not affected by intermittency depends on the 
first power of e) p2| . 

In Figure we plot the first moments of inverse dou- 
bling time (In) compensated with the Kolmogorov scaling 
R~ 2p / 3 . The quality of the scaling is remarkable, espe- 
cially if compared with the standard statistics of Figure 0. 
This allows to detect small deviations from dimensional 
scaling. Indeed, a closer inspection of Figure H reveals 
that the compensation is not perfect, the deviation be- 
ing more evident for higher moments; this indicates the 
existence of Lagrangian intermittency. 

Figure |5j shows some moments of the inverse doubling 
time, now compensated with best fit exponents /3 p . The 
improvement with respect to Figure |4] demonstrates that 
the exponents in (0) are corrected in comparison to di- 
mensional prediction. From the doubling time exponent 
/3„ we can obtain the Eulerian exponents C, by inverting 
(g) . The result shown in the inset of Figure [| gives a set of 
exponents Cp which are remarkably close to "standard" 
structure function exponents in fully developed turbu- 
lence. We stress that, at the present resolution, the scal- 
ing of the Eulerian structure function is rather poor, thus 
a precise determination of C P is possible only using indi- 
rect analysis, such as the ESS technique p3fl . 

Let us summarize our findings. We have performed 
direct numerical simulations of a three-dimensional tur- 
bulent flow and concentrated on the problem of parti- 
cles' dispersion. The overall dispersion behavior is well- 
described by the Richardson's equation, although some 
deviations (mostly caused by the finite-Reynolds nature 



of the simulations) are evident. The use of fixed-scale 
statistics (doubling-time distribution) instead of fixed- 
time ones removes to a large extent these restrictions, 
and gives a possibility to evaluation the Richardson's con- 
stant very accurately. Its value is Ci ~ 0.55, in a close 
agreement with recent experimental findings. The dis- 
cussion of the inverse moments of the doubling-time dis- 
tributions unveils the role of Lagrangian intermittency in 
the two-particle dispersion. The values of the Lagrangian 
scaling exponents are connected with the Eulerian struc- 
ture function exponents Q p . The values of Cp obtained 
from the separation statistics are remarkably close to 
standard ones, known for fully developed turbulence. In 
the next future it will be probably possible to have exper- 
imental Lagrangian trajectories in high Reynolds number 
flows [£4| . It would be extremely interesting to check our 
findings in real fluid turbulence. 
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FIG. 3. Probability distribution function of relative sepa- 
rations at three different times. The continuous line is the 
Richardson prediction (fa), the dashed line is the Gaussian 
distribution proposed by Batchelor. 
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FIG. 1. Average energy spectrum E(k). The dashed line 
has the Kolmogorov slope —5/3. In the inset it is shown the 
compensated third order longitudinal structure function. The 
dashed line represents the "4/5" law. 



FIG. 4. First moments of the inverse doubling time 
((1/T(R)) P ) compensated with Kolmogorov scaling R~ 2p ' 3 . 
Deviations from dimensional compensation are evident, in 
particular for p — 4. In the inset we plot the compensated 
mean doubling time according to (pi) together with the esti- 
mate corresponding to C2 — 0.55. 
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FIG. 2. Relative dispersion R 2 (t) versus time t. The 
dashed line is the Richardson i 3 law. In the inset we show the 
compensated plot R (t)/(et ) which should give the Richard- 
son constant Ci. Because of the strong oscillation, a precise 
estimation of C2 is difficult. 
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FIG. 5. First moments of the inverse doubling time 
{(1/T(R)) P ) compensated with best fit exponent (5 P . Observe 
the improvement in the compensation with respect to Fig- 
ure W. In the inset we plot the structure function exponents 
estimated from £ p = p + /3 P . The line is the She-Leveque 
parameterization. 
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